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ABSTRACT 

PHOEBE (PHysics Of Eclipsing BinariEs) is a modeling package for eclipsing bi¬ 
nary stars, bnilt on top of the widely used WD program (Wilson & Devinney 
1971). This introductory paper overviews most important scientihc extensions 
(incorporating observational spectra of eclipsing binaries into the solution-seeking 
process, extracting individual temperatures from observed color indices, main- 
sequence constraining and proper treatment of the reddening), numerical inno¬ 
vations (suggested improvements to WD’s Differential Corrections method, the 
new Nelder & Mead’s downhill Simplex method) and technical aspects (back-end 
scripter structure, graphical user interface). While PHOEBE retains 100% WD 
compatibility, its add-ons are a powerful way to enhance WD by encompassing 
even more physics and solution reliability. The operability of all these extensions 
is demonstrated on a synthetic main-sequence test binary; applications to real 
data will be published in follow-up papers. PHOEBE is released under the GNU 
General Public License, which guarranties it to be free, open to anyone interested 
to join in on future development. 

Subject headings: methods: data analysis — methods: numerical — binaries: 
eclipsing — stars: fundamental parameters 
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1. Introduction 

With the ever-growing computer power, numerical models built to analyze acquired 
eclipsing binary data are gaining both on accuracy and complexity. The motivation is clear: 
due to their unique geometrical and kinematic properties, eclipsing binaries (EBs) give full 
physical insight into their structure, distance and evolution stage of their coeval components. 
In the last 40 years the EB held was overwhelmed by many approaches to solution seeking; 
Kallrath & Milone (1999) give an overview of most important ones. The widely used WD 
code (Wilson & Devinney 1971) underwent many expansions, improvements and hue-tuning 
(Wilson & Soha (1976); Wilson (1979, 1990); Milone et ah (1992); Kallrath et al. (1998); 
Van Hamme & Wilson (2003) and many others), which hrmly established it as the most 
prominent software available for EBs. 

So why would one build yet another modeling program? The answer is simple: one would 
not. Tackling same old problems all over again does not make sense; rather, one builds on 
basis of what has already been done. This is what our ehort is all about: to create a modeling 
package built on top of the Wilson-Devinney code, introducing new enhancements to where 
WD was dehcient, while still pertaining 100% WD compatibility. Enhancements include 
new physics (proper handling of color indices and therefore temperatures in absolute units, 
interstellar reddening ehects), existing minimization scheme add-ons (stability and conver¬ 
gence improvements) and new minimization schemes aiming to fully automate hrst steps of 
solution-seeking (an issue of utmost importance for ambitious space scanning missions like 
Gaia (Perryman et al. 2001)). We discuss main characteristics of this new package called 
PHOEBE: PHysics Of Eclipsing BinariEs. 

This paper introduces the formalism PHOEBE is built on and demonstrates its capabilities 
on synthetic binary data that are described in Section 2. Sections 3 and 4 give a detailed 
overview of computational and physical extensions. Section 5 discusses further work to be 
done and explains future vision of this project. Technical details on PHOEBE availability and 
license, back-end logic and structure, and front-end interface are given in the Appendix. 


2. Building a test binary star 

To demonstrate innovations PHOEBE brings to the EB held, a synthetic binary model 
is created. Testing the methods against a synthetic model may seem artihcial, but the 
obvious advantage of knowing the right solutions is the only true way of both qualitative 
and quantitative assessment. Some preliminary results of using PHOEBE on true observations 
were already presented by Prsa (2003). Full-hedged demonstration of PHOEBE capabilities 
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both for individual stars and large data-sets will be published in series of follow-ups to this 
paper shortly. 

Our synthetic binary consists of two main-sequence F8 V-Gl V components with their 
most important orbital and physical parameters listed in Table 1. It is a partially eclipsing 
detached binary with only slight shape distortion of both components (i?i,poie/.Ri, point = 
0.974, i?2,poie/7?2,point = 0.979). Light curves are generated for Johnson B and V passbands 
in 300 phase points with Poissonian scatter ranging from ay = 0.005 to ay = 0.025 at 
quarter-phase magnitude my = 10.0. Radial velocity (RV) curves are generated in 50 phase 
points with Gaussian scatters ranging from ctrv = Ikms”^ to ctrv = 25 km s~^. Light 
curves in B and V with ay = 0.015 and both RV curves with (Trv = 15km s“^ are depicted 
in Fig. 1. 

This model binary will be used for demonstrating all PHOEBE’s capabilities that are novel 
to the field of EBs. 


3. Solving the inverse problem for eclipsing binaries 

The underlying WD code is composed of two parts: the LC program for computing light 
and RV curves and the DC program for solving the inverse problem (Wilson 1993). PHOEBE 
introduces several optimizations to the DC method and adds to generality by implementing 
a new minimization method: Nelder & Mead’s downhill Simplex. 


3.1. Suggested optimizations to WD solving method 

WD’s DC code, as the name suggests, uses Differential Gorrections (DG) method comple¬ 
mented by the Levenberg-Marquardt algorithm to solve the inverse problem (Wilson 1993). 
It is especially suited for EBs and is one of the fastest codes around. In cases when the 
method does not converge, the Method of Multiple Subsets (MMS) may be used to relax the 
system to the nearest minimum (Wilson & Biermann 1976). 

A DC program reads in a user-supplied input hie consisting of a) a set of initial param¬ 
eters that dehne physical and geometrical properties, b) observational data and c) switches 
that dehne the way a minimization algorithm is run (refer to the booklet by Wilson & Van 
Hamme (2003) accompanying WD code for details on DC input hies). Within one iteration, 
the values of parameters set for adjustment are improved and returned for user inspection. 
In case of convergence, the user manually resubmits the new parameter set to the next it¬ 
eration. The measure of the quality of the ht (the cost function) is the sum of squares of 
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Table 1. Physical parameters of the F8 V-Gl V test binary star. Spectral type - 
temperature relation taken from Lang et ah (1992). 


Parameter [units] Binary 

F8 V G1 V 


Po [days] 


1.000 


a [Rq] 


5.524 


q = m 2 /mi 


0.831 


^ r] 


85.000 


[km s ^] 

Teff [K] 

6200 

15.000 

5860 

L [Lo] 

2.100 


1.100 

M [Mq] 

1.239 


1.030 

R [Rq] 

1.260 


1.020 

G H- 

5.244 


5.599 

log {g/go) [-]‘’ 

4.33 


4.43 

Xb [H'" 

0.818 


0.833 

VB [-]^ 

0.203 


0.158 

Xv [-]^ 

0.730 


0.753 

yv [-]" 

0.264 


0.242 


^Unitless effective potentials dehned by Wilson (1979). 

= 1cm s“^ is introduced so that the logarithm acts on a dimensionless variable. 

^Linear (x) and non-linear (y) coefficients of the logarithmic limb darkening law for Johnson 
B and V passbands, taken from Van Hamme (1993). 
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Fig. 1.— F8 V-Gl V test star data. Light curves are computed for Johnson B (Filed 
dots) and V (empty dots) passbands in 300 phase points with ay = 0.015 (left panel). RV 
curves are computed in 50 phase points with o-rv = 15km s“^ (upper right panel); eclipse 
proximity effects are turned off. Star plot is computed at quarter phase, cross denotes the 
center-of-mass (lower right panel). 
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weighted O — C residuals. 

WD’s list of more than 30 adjustable parameters includes passband luminosities L\ for 
i light curves (with their WD name HLA), that have a unique property of linearly scaling 
the level of light curves. DC (or any other minimization algorithm) fits these luminosities the 
same way it hts all other physical parameters: softly. This means that within one iteration, 
the values of L\ are not fully adjusted, only improved. Since Lj’s determine vertical offset 
of light curves, this raises two specific problems: 1) The soft change of L\ in every iteration 
step causes artificial changes of other physical parameters: rather than fitting the shape 
of the data curve, other parameters £t the discrepancy between the model and the data, 
induced by the softness of L\ fit. It is like driving a very old car on a very bnmpy road - 
each bnmp on the road causes wobbling of the whole car with slow attenuation. 2) Changes 
of adjusted parameters calculated by DC will properly contribute to the cost function only if 
the model is aligned with the data: the average O — C valne must be approximately 0. This 
alignment is governed by L\ for light curves. If this alignment is not computed correctly, the 
cost function is misleading DC instead of aiding it. This causes under-estimation of formal 
errors dne to L\ softness error propagation and even convergence problems. 

PHOEBE solves this problem by snpplying an option to compute Lj’s instead of minimiz¬ 
ing them, thus increasing their stiffness with respect to other parameters. The alignment is 
calculated so that the average O — C value is exactly 0. The time cost of this computation 
is not only negligible, it actually speeds up the overall algorithm, since the dimension of the 
parameter subspace submitted to DC is reduced. Fig. 2 demonstrates the iteration seqnence 
with the original method (left) and the proposed method (right) for a case of 7 simnlta- 
neously fitted parameters displaced by at most 50% from their true value. In the latter 
case parameters converge quickly and in a smooth fashion. Similar simulations that test 
convergence behavior in cases when both temperatures are fitted or when other individual 
parameters are kept constant have also been performed; they accord or even amplify the 
conclusion of Fig. 2 and their results are thus omitted on account of brevity. Note however 
that stiffening Lj’s does not guarantee convergence to the global minimum, it only solves 
the inverse problem more efficiently. It should also be stressed that calculating Lj’s instead 
of fitting them might not always affect convergence as noticeably, particularly in cases where 
relative corrections of parameter valnes are small. 

By calculating Lfs instead of fitting them, the criterion is not used and the corre¬ 
sponding formal errors of Lj’s are not calculated. To obtain them, one would simply revert 
from calculating to fitting Lj’s at the very end of the minimization process and submit them 
to the final iteration of the DC. 
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Fig. 2.— Soft vs. stiff curve leveling. Iteration sequence for 7 physical parameters of our test 
binary for soft L\ scheme (left) and stiff L\ scheme (right). The x-axis is given in log-scale 
to amplify the part where the impact of stiffening is largest. Unity on g-axis corresponds to 
parameter’s initial value. Li and L 2 are passband luminosities in B and V hlter, respectively, 
i is the system inclination, L2i and L 22 are gravity potentials, q is the mass ratio and Ti and 
T 2 are surface temperatures. Temperature Ti is kept constant throughout the £t, simulating 
the usual practice of determining one temperature and htting the other. 


Systemic velocity v^. The levels of RV curves are determined by the systemic velocity 
v^\ changing it vertically shifts those curves. Although is not as correlated with other 
parameters as is the case for L|’s and the problem is thus not as severe, alignment between 
the model and the data is still crucial. PHOEBE allows calculation following the same logic 
as before for L^’s - by demanding that the average O — C value is exactly 0. 


Limb darkening coefficients. The native WD code supports linear, logarithmic and 
square root limb darkening (LD) laws. Their coefficients primarily depend on the given 
passband, effective temperature, gravity acceleration log {g/go) and metallicity [M/H]. WD 
does not constrain the choice of these coefficients, so people have traditionally used LD tables 
computed by, e.g.. Van Hamme (1993) or Claret (2000). 

Following a similar argument to the one mentioned before for L*^’s and PHOEBE im¬ 
plements an optional dynamical LD computation. After each iteration that induces changes 
to any of the Teg, log {g/go), [M/H] or related parameters, the LD coefficients need to be 
modified accordingly. PHOEBE uses Van Hamme (1993) tables for this purpose, dynamically 
reading out tabulated values and linearly interpolating to obtain proper values automati¬ 
cally. The implications are not as severe as for the L^’s and because LD contributions 













are orders of magnitude smaller and insensitive to small changes in the above mentioned 
parameters. 


3.2. New minimization algorithms 

The main driving force of any binary minimization algorithm is its ability to solve the 
inverse problem as accurately and as quickly as possible. WD’s DC algorithm is very fast 
and works well if the discrepancy between the observed and computed curves is relatively 
small, but it can diverge or give physically implausible results if the discrepancy is large. 
While this dehciency is usually not a severe problem when analysing individual EBs (one 
can always obtain a reasonable set of starting parameters by calculating a few initial light 
and RV curves), its impact when dealing with huge data-sets (such as hundreds of thousands 
of light curves that will be obtained by Gaia) may be a blocker. To overcome this, and to 
assist in initial steps of solution-seeking, a complementary minimization scheme to DC is 
proposed. 


Nelder &; Mead’s downhill Simplex. Two main dehciencies of DC are especially strik¬ 
ing. 1) The main source of divergence and the loss of accuracy in DC is the computation of 
numerical derivatives of the cost function with respect to parameters set for adjustment. 2) 
Once DC converges, there is no ready way of telling whether the minimum is local or global; 
the method cannot escape. The latter problem affects most minimization algorithms that 
have been applied to EBs. 

To circumvent these two problems, PHOEBE implements Nelder & Mead’s downhill Sim¬ 
plex^ method (Nelder & Mead 1965), hereafter NMS. Since NMS does not compute deriva¬ 
tives but relies only on function evaluations, it cannot diverge. The basic form of NMS 
applied to a WD implementation was first proposed by Kallrath & Linnell (1987). PHOEBE 
goes a step further and adapts the method specihcally to EBs. Eirst tests of PHOEBE’s NMS 
implementation on photometric data that are expected to be obtained by Gaia (Prsa & 
Zwitter 2005) are very promising. 

NMS acts in n-dimensional parameter hyperspace. It constructs n vectors pj from the 
vector of initial parameter values x and the vector of step-sizes s as follows: 

Pi (^0) ^1) ■ ■ ■ t T Sj, . . . , Xfi') (1) 


^Nelder & Mead’s downhill Simplex should not be confused with linear or non-linear programming algo¬ 
rithms, which are also referred to as Simplex methods (e.g. Press et al. 1992). 
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These vectors form (n + 1) vertices of an n-dimensional simplex. During each iteration 
the algorithm tries to improve parameter vectors p* by modifying the vertex with the high¬ 
est function value by simple geometrical transformations: reflection, reflection followed by 
expansion, contraction and multiple contraction (Galassi et ah 2003). Using these transfor¬ 
mations, the simplex moves through parameter space towards the closest minimum, where 
it contracts itself. Fig. 3 shows the number of iterations required for the NMS to converge 
from different starting points within 10“^ fractional accuracy. The data are those of our test 
binary given in Fig. 1. 

This basic form of NMS is unconstrained, which means that parameters may assume any 
value regardless of their physical feasibility. The NMS implemented by PHOEBE optionally 
enables semi-constrained or fully constrained minimization by imposing limits to several or all 
adjusted parameter values. Additionally, heuristic scan, parameter kicking and conditional 
constraining enable NMS to efficiently escape from local minima. 


3.3. Heuristic Scan 

EB minimization algorithms, including even NMS with its property of guarranteed con¬ 
vergence, can be stuck in a local minimum, particularly since parameter hyperspace in vicin¬ 
ity of the global minimum is typically very flat, with lots of local minima. In addition, global 
minimum may be shadowed by data noise and degeneracy. 

Heuristic scan is an enhancement method to any minimization algorithm (DC, NMS, 
...) that selects a set of starting points in parameter hyperspace and starts the minimiza¬ 
tion from each such point. The user defines how starting points are selected - they may 
be gridded, stochastically dispersed, distributed according to some probability distribution 
function (PDF) etc. The algorithm then sorts all solutions by the cost function (the 
for example) and weights the obtained parameter values accordingly: heuristic runs with 
smallest values of the cost function correspond to the deepest minima and should thus be 
most weighted - they are most suitable candidates for the global minimum. 

The weighted values of adjusted parameters are then put into histograms, from which the 
mean and standard deviation of parameter values are calculated. These estimates are truly 
statistical, since they do not depend on formal errors of the numerical method. Fig. 4 shows 
an example of such histograms for the effective temperature ratio r = T 2 /T 1 . Heuristic scan 
results for this particular example are virtually insensitive to observational data accuracy: 
for three signihcantly different cases (labelled best, medium and worst quality data on Fig. 4), 
the outcome of the histogram £t is approximately the same. Histograms for other parameters 
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have somewhat larger standard deviations because of degeneracy: obtained inclination for 
medium quality data is 85.6° ± 1.07° (compared to the true value i = 85°) and gravitational 
potential fli is 5.44 ± 0.27 (compared to the true value fli = 5.244). It should be noted 
that reliable statistics implies many starting points of heuristic scan, which in turn implies 
signihcant prolongation of the algorithm computation time: each additional scan linearly 
contributes to the time cost. 

Because of data noise and degeneracy, the global minimum is essentially never a single 
point (with its corresponding uncertainty), it is actually a region (with its corresponding 
uncertainty) in parameter hyperspace. Such a region encompasses many adjacent minima, 
the depths of which are physically indistinguishable - a single observed data point with its 
individual weight may change the identity of the deepest minimum within that region. To 
identify these regions, PHOEBE computes convergence tracers - selected 2D cross-sections of 
the parameter hyperspace, tracing parameter values from each starting point, iteration after 
iteration, all the way to the converged solution. Attractors - regions that attract most con¬ 
vergence traces - within these cross-sections reveal parameter correlations and degeneracy. 
Inspecting such convergence tracers offers additional insight on the quality and integrity of 
the solution. Local minima in context of convergence tracers are those that lay outside of 
the deepest attractor(s); those are the ones that need to be identihed and escaped from. 

A particularly troublesome degeneracy is the one between the inclination and either of 
the effective potentials of the two stellar components (which act on behalf of components’ 
radii). Fig. 6(a) shows the i-Qi convergence tracer computed for our test binary. The 
correlation between i and Di is evidently very flat at i ~ 85°, which may be easily understood: 
the model is able to compensate smaller inclinations by enlarging the radius of the star and 
vice versa. Therefore we should not trust light curve analysis to disentangle these parameters 
by itself - additional constraints are needed. This issue will be further discussed in Section 
4. 


3.4. Parameter Kicking 

Another possible approach to detect and escape from local minima is to use a stochastic 
method such as Simulated Annealing (SA). However, such methods are notoriously slow. 
Thus, instead of full-featured SA scan, a simple new procedure has been developed that 
achieves the same effect as stochastic methods, but in signihcantly shorter time. The idea is 
as follows: whenever a minimum is reached within a given fractional accuracy, the algorithm 
runs a globality assessment on that minimum. If we presume that standard deviations ak 
of observations are estimated properly and that they apply to all data points, we may use 
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Fig. 3.— Histogram of the number of iterations required for NMS convergence within 10“^ 
fractional accuracy. The Probability distribution function (PDF) exhibits a maximum at 
~75 iterations. Due to extremely fast convergence in the first few steps, the number of 
iterations is in practice insensitive to the selection of the initial starting point in parameter 
hyperspace; the required number of iterations is dominated by convergence behavior in the 
” Minima valley”. 



Fig. 4.— Temperature ratio histogram obtained as a result of heuristic scan. Plots show 
T = T 2 /T 1 PDFs for three different observational datasets: ctlc = 0.005, (Try = 5kms“^ 
(left), (Tlc = 0.015, (Try = 15km s“^ (middle) and ctlc = 0.025, ctry = 25km s“^ (right). 
First and last bins hold all other outlying points. Heuristic scan is practically insensitive 
to the observational data accuracy as long as there are sufficient data points to determine 
both eclipse depths. Obtained values of temperature ratios are purely statistical and may 
be compared to the true value of r = 0.9452. 
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them for weighting: 

M 1 ^ 

^k = '^ - Vif = — X] Wi{Xi - Vif, 

1=1 ^ t=l 


( 2 ) 


where index i runs over M measurements within a single data-set and index k runs over N 
data-sets (photometric and RV curves); Xi are observed data points, yi are calculated data 
points and tCj’s are individual weights. Since the weighted variance is given by: 


s 


2 

k 



^Wi{Xi - Vif, 


(3) 


we may readily express x\ as: 


and the overall x^ value as: 


xl 






(4) 

(5) 


If (Jk are realistic, the ratio Sk/ck is of the order unity and of the order iVtot = ^k- 

This we use for parameterizing values: 


A := {x^/Ntot) ■ quantization. 


( 6 ) 


Parameter kicking is a way of knocking the obtained parameter-set out of the minimum: 
using the Gaussian PDF, the method randomly picks an offset for each parameter. The 
strength of the kick is determined by the Gaussian dispersion Ukick, which depends on the 
minimum globality assessment parameter A. If A is high, then the kick should be strong, but 
if it is low, i.e. around A ~ 1, then only subtle perturbations should be allowed. Experience 
shows that a simple expression such as: 


_ 0.5A 

t^kick - ^ 


(7) 


works very efficiently in case of partial eclipses. This causes Ukick to assume a value of 0.5 for 
lOa offsets and 0.005 for la offsets, being linear in between. Note that this dkick is relative, 


i.e. given by: 


cr, 


abs 

kick 


X(J 


rel 

kick’ 


( 8 ) 


where x is the value of the given parameter. When convergence within the given fractional 
accuracy is reached, parameters are kicked with respect to the depth of the minimum and 
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the minimization is restarted from displaced points. The influence of consecutive parameter 
kicking with NMS is depicted in Fig. 5; it is shown that out of all heuristic scans only ~30% 
initially converge to within 1% of optimal value of A, whereas this percentage steadily grows 
to ~60% after three kicks. Figs. 6(b)-(d) show signihcant improvement to the solution 
introduced by these consecutive kicks. Parameter kicking is able to quickly escape from 
local minima and thus rapidly increase convergence efficiency of the whole NMS method. A 
down-side of parameter kicking is the time cost: each additional kick linearly adds to the 
overall execution time. For thorough discussion and details on benchmarking please refer to 
PHOEBE accompanying documentation. 

The idea behind the NMS implementation is not to replace DC, but rather to com¬ 
plement it. DC is created for interactive usage and converges in discrete steps that need 
monitoring. NMS on the other hand aims to automate this process so that intermediate 
monitoring is no longer necessary. DC is one of the fastest methods (WD’s DC in particular, 
since it is optimized for EBs), but may easily diverge. At the expense of speed, NMS is one of 
the most robust algorithms for solving non-linear minimization problems and never diverges. 
Finally, both DC and NMS suffer from degeneracy and may become stuck in local minima. 
To overcome this, both methods are complemented by heuristic scan and parameter kicking. 
These differences in intent make a combination of the two methods a powerful engine for 
solving the inverse problem. 


4. Extended set of physical constraints 

WD’s extensive list of more than 30 adjustable parameters is an overwhelming indica¬ 
tor of how sophisticated the model has become in 35 years of development. Nevertheless, 
accuracy is crucial for a model to describe such a wide diversity of intrinsically different 
binaries. An accurate model should contain all relevant physical contributions for which the 
governing laws are well-known. We start the discussion by introducing new physical ties and 
constraints to parameter extracting schemes that are implemented in PHOEBE. It should be 
stressed that all these constraints are optional and it is up to the user to select the ones that 
are of interest. 


4.1. Color indices as indicators of individnal temperatures 

One of the main difficulties of modeling EBs is accurate determination of individual 
temperatures of both components. Frequent practice in literature is to assume the temper- 



ature of one star (e.g. from spectra or color indices) and fit the temperature of the other 
star. This approach is often inadequate, particularly for binaries with similar component 
temperatures and luminosities: in such cases, the contribution of both components to the 
system luminosity is signihcant and it is difficult to accurately estimate the contribution of 
only one star in advance. 

Before we propose a method capable of providing individual temperatures from standard 
photometry observations without any a priori assumptions, it proves useful to introduce a 
concept of effective temperature of the binary. A binary may be regarded as a point-source, 
the effective temperature of which varies in time. To this effective temperature contribute 
both components according to their sizes and individual temperatures, and the inclination. 
Effective temperature of the binary is directly revealed by the color index, so its observational 
behavior is well known. If a model is to accurately reproduce observations, the composite of 
contributions of both components must match this behavior. 

The observational light curve quantity (dependent variable) WD works with is flux, 
scaled to an arbitrary level (which could also be in absolute physical units, i.e. W/m^ per 
wavelength interval). The model adapts to this level by determining the corresponding 
passband luminosity L\, one for each passband. However, these passband luminosities are 
completely decoupled from one another, so any color information that might have been 
present in the data is discarded. Since the effective temperature of a binary is observationally 
revealed by its H H (or any other suitable) color index^, some of the relevant temperature 
information is lost. Transformation to fluxes in absolute units would not suffice for properly 
determining the corresponding passband luminosities - one needs a physical relation between 
those luminosities. Neglecting this additional relation may result in discordant colors between 
the temperatures obtained by the £t (assuming that Ti is a priori known) and the ones 
determined by the binary’s effective temperature. This relation is nothing else than the 
color index and may thus be accurately determined from observations. 

In the last decade substantial effort was made to scan the sky for standard stars to 
be used for photometric calibration: Landolt (1992) covering celestial equator, Henden & 
Honeycutt (1997) and Bryja & Sandtorf (1999) covering fields around cataclysmic variables, 
Henden & Munari (2000, 2001) covering fields around symbiotic binaries, to name just a 
few. These efforts help overcome the problem of small CCD helds with respect to all-sky 
photometry, since in many fields there are now cataloged standard stars that may be used 
to extract color indices for EBs. In context of PHOEBE, this means that using measured 
color indices as additional information is plausible even if the data were not obtained under 


^Useful relations among color indices are given in Caldwell et al. (1993). 
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photometric conditions. 

PHOEBE initially regards L*^’s as simple level-setting quantities - physical context comes 
in only after the color index constraint is set. For the sake of simplicity, consider that input 
observational data are supplied in magnitudes (such that the color indices are meaning¬ 
ful) rather than fluxes. PHOEBE input data in individual passbands should not be scaled 
arbitrarily; that is PHOEBE’s job. 

The native type that PHOEBE works with is inherited from WD, which is flux. PHOEBE 
uses a single, passband-independent parameter mo to transform all light curves from magni¬ 
tudes to fluxes. The value of mo is chosen so that the fluxes of the dimmest light curve are 
of the order of unity. It is a single quantity for all light curves, which immediately implies 
that the magnitude difference, now the flux ratio, is preserved; hence, the color index is 
preserved. If the distance to the binary is known (e.g. from astrometry), mo immediately 
yields observed luminosities of the binary; intrinsic luminosities are obtained if the color 
excess E{B — V) is also known. 

This is where physics comes in; from such set of observations, the calculated Lj’s are 
indeed passband luminosities, the ratios of which are the constraints we need: passband 
luminosities of light curves are now connected by the corresponding color indices. Once the 
color constraint is set, PHOEBE makes sure that the ratio between Lj’s is kept constant. 

Now that the color indices are preserved, effective temperature of the binary may be 
obtained from a color-temperature calibration. PHOEBE uses updated Flower (1996) tables 
with coefficients given in Table 2. It should be stressed that the color constraint is applicable 
only if the data are acquired on (or properly transformed to) a standard photometric system. 

Applying the color constraint, effective temperatures of individual components may be 
readily disentangled by the minimization method. The method is now able to hud only those 
combinations of parameters that preserve effective temperature of the binary and hence the 
color index. Since the relation between effective temperatures of individual components is 
fully determined by the light curve shape (dominantly by the primary-to-secondary eclipse 
depth ratio) and since the sum of both components’ contributions must match the effec¬ 
tive temperature of the binary, the color-constrained minimization method yields effective 
temperatures of individual components without any a priori presumptions. 

Let us demonstrate this concept on our test binary. Calculating passband luminosities 
from medium quality observations (ctlc = 0.015, ctrv = 15 km s“^) yields Lb/L y = 0.592 ± 
0.006. Transforming this into magnitudes yields the color index B — V = 0.57 ±0.01, which 
in turn yields effective temperature of the binary to be = 6 002 K ± 40 K. The relation 
between both individual temperatures from the ratio of eclipse depths is well determined 
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A 


Fig. 5.— A-histogram for initial heuristic scan and three consecutive parameter kicks. The 
success of parameter kicking is obvious, since after only three consecutive kicks the percentage 
of scans that converge within 1% of the optimal value of A (in case of proper a^s A = 
1) is doubled from ~30% to ~60%. As is shown by Prsa & Zwitter (2005), parameter 
kicking proves to be even more efficient in case of exclusively photometric observations, with 
improvement from ~15% to ~75% in convergence after three consecutive kicks. 


Table 2. Coefficients of the empirical Teff {B — V) relation given by the 7th degree 
polynomial £t ~ (Flower, private communication). The second 

column applies to main-sequence stars, sub-giants and giants, the third column applies to 

supergiants. 


Coefficient; 

V, IV, III, II 

1 

C'o 

3.979145 

4.012560 

C*! 

-0.654992 

-1.055043 

^2 

1.740690 

2.133395 

C's 

-4.608815 

-2.459770 

C'4 

6.792600 

1.349424 


-5.396910 

-0.283943 


2.192970 

— 

C'7 

-0.359496 

— 
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Fig. 6.— Convergence tracer for cross-section. This particular case is a notorious 
example of very difficult-to-handle correlation between the inclination and effective potentials 
(hence the radii) of both components (only correlation is depicted for brevity). Individual 
plots denoted with letters (a) through (d) show the result of NMS heuristic scan from zero 
to three consecutive parameter kicks. Cross-hairs mark the position of the true minimum. 
Attractors are symmetric to i = 90°, but still very flat at i ~ 85° to 90° interval, which 
means that the obtained NMS solution should not be blindly trusted; rather, additional 
constraining is needed. 
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(c.f. Fig. 4 yielding T 2 /T 1 = 0.946±0.003), disentangling effective temperatures of individual 
stars to be Ti = 6 190 K ± 52 K and T 2 = 5 880 K ± 51 K. Comparing these values to true 
values Ti = 6 200 K and T 2 = 5 860 K is very encouraging. 

4.2. Spectral energy distribution as independent data source 

Traditionally, spectral energy distributions (SED) have been used only indirectly, e.g. 
for extracting radial velocities or determining effective temperatures. Several recent studies 
of individual EBs have shown that including flattened SEDs may be used as external check 
of the model solution (see Siviero et ah (2004) and Marrese et ah (2004) for examples), 
where individual spectral lines of Echelle spectra are compared with Kurucz (1998) model 
atmospheres. 

Since the Kurucz’s model atmosphere program runs only under VAX/VMS in its dis¬ 
tributed form, several databases of precomputed spectra have been assembled for practical 
use (e.g. Zwitter et ah (2004) covering the spectral range 765-875nm, Murphy & Meiksin 
(2004) covering 300-1000nm, Munari et ah (2005) covering 250-1050nm). Such databases 
bring stellar atmospheres to non-VAX/VMS equipped users. In recent years signihcant effort 
has been made to port Kurucz’s model atmospheres code to Linux (see e.g. CCP7 initiative 
at http://www.stsci.edu/software/CCP7, Sbordone et al. (2004) and others). Such initia¬ 
tives enable users to include SED data in solving the inverse EB problem. PHOEBE already 
takes a step in that direction by using a synthetic spectra database to test whether flattened, 
wavelength-calibrated spectra match synthetic spectra within a given level of significance. 

One very important caveat that should be stressed: it is not feasible to compare obser¬ 
vational SEDs to synthetic SEDs over the full spectral range. The problems occur because of 
Earth’s atmosphere (signihcant parts of the spectrum are dominated by telluric lines, which 
the model does not handle). By default, PHOEBE uses the Zwitter et al. (2004) grid of 61196 
synthetic spectra covering the 765-875 nm interval at a resolving power i? = 20 000. A 
simple interpolation may be used to obtain the spectrum characterized by any combination 
of Teff, log (s'/s'o), [M/H] and Urot with the accuracy better than 25 K in temperature, 0.05 
dex in log {g/go) and metallicity and 1 km s~^ in rotational velocity. These uncertainties 
are smaller than the uncertanties of the Kurucz’s model for parameters of our test binary, 
so interpolation does not induce any systematic errors. 

To demonstrate current level of SED implementation in PHOEBE, consider again our test 
binary. Parametric vectors (Teg, log {g/go), nrot)i ,2 of both EB components are determined by 
the model solution from photometric and RV data. These are used to obtain synthetic spectra 



by linear interpolation in Te^, log {gjgo) and Wrot from the grid. For the ’’true” simulated 
spectrum, solar abundances ([M/H] = 0.0), corotation (uroti = 64km s“\ = 52km s“^) 

and microturbulence Vturb = 2km s“^ are assumed. Effective spectrum of the binary is 
computed in out-of-eclipse phase (e.g. quarter-phase) by Doppler-shifting and convolving 
the spectra of the two stars. Fig. 7 shows an example of a quarter-phase spectrum of the 
test binary. 

After each solution of the NMS heuristic scan, a synthetic spectrum is built^ from that 
solution. It is then compared with the ’’true” spectrum by the cost function. Effective 
temperature is the dominant parameter that governs SED shape, but this is of little use for 
our case: recall from Fig. 4 and the discussion on color indices that individual temperatures 
are well determined from photometry alone. Rather, our solution suffers from degeneracy in 
effective potentials Di, 142 and inclination i (Fig. 6). It would be benehcial if the SEDs could 
break this degeneracy. Since the mass ratio and semi-major axis of the model are effectively 
held constant by the RVs, Di and 142 depend only on the radii of individual components. 
Thus, different D’s imply different log {g/go) and, by assuming corrotation, also v^ot- Fig. 8 
shows the Vroti-'yrot 2 cross-sectiou, demonstrating that, as we hoped, the SED analysis indeed 
constrains the solution to smaller intervals for Wj-oti and Wrot 2 ; thus smaller intervals for Di 
and D 2 - 

The Wroti-nrot 2 cross-sectiou may sometimes do even more than only break the degen¬ 
eracy between Di and 122 - If the radii are well determined, e.g. by total eclipse geometry, 
such analysis yields synchronicity parameters F\ and F 2 , since the only way to compensate 
the change in rotational velocities for any predetermined radii is to break the corotation 
presumption. This may be especially important in analysis of well detached systems, as 
demonstrated by Siviero et al. (2004). 

It should be noted that there is no support for extracting Teg, log {g/go), [M/H] or 
Urot from spectra at the moment, only a weighted f^st is done to conhrm or reject the 
particular set. As such, the current implementation forms the base of spectral analysis 
for EBs, but it still does not contribute fully to minimization. Once we are capable of 
building stellar spectra without presuming spherically-symmetrical stars in LTE, full SED 
will be introduced to the minimization process as well. Such a scheme will have to weight 
properly individual wavelengths, since there is much less information in the continuum of 
the spectrum than it is, in example, in central parts and wings of spectral lines. However, 
even the present implementation of SED analysis hnds the values of physical parameters 
which have not usually been attainable by light and RV curve analysis, namely metallicity 


•^At present the spectrum may be generated for any orbital phase outside of eclipses. 
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Fig. 7.— Synthetic spectrum of the test binary at quarter phase. The spectrum is built by 
Doppler-shifting and convolving individual component spectra obtained by linear interpola¬ 
tion in Teff and log {g/go) from precomputed stellar spectra tables by Munari et al. (2005). 
The inset magnihes a part of the spectrum corresponding to the Gaia RVS wavelength range, 
which is covered by the Zwitter et al. (2004) database. The strongest lines in the inset are 
split (revealing the binary nature of the object) and are due to Ca II (849.80 nm, 854.21 nm 
and 866.94 nm). 
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Fig. 8.— Result of the comparison of the ’’true” spectrum against Zwitter et ah (2004) 
database. Out of all cross-sections, Uroti-'yrot 2 cross-section is most interesting, because it 
helps break the degeneracy between effective gravitational potentials Oi and 02- The levels 
of gray in the mesh are linear in and denote the quality of the £t: black color corresponds 
to the best fit, white color corresponds to the worst fit. Cross-hairs denote the position of 
the true values of rotational velocities. 



and rotational velocity (see Terrell et al. 2003). 


4.3. Main sequence constraints 

In cases where SED observations are not available, or where they are used only to 
extract RVs, the degeneracy among parameters may still be so strong that neither heuristic 
scan nor parameter kicking can break it. In such cases we stand no chance of obtaining any 
satisfactory solution without further constraining the modeled binary. 

WD features 8 modes of operation that determine the morphology of the binary. By 
deciding on the mode of operation, the user imposes a set of physical constraints; for example, 
both components of over-contact systems have equal potentials. PHOEBE refers to these 
constraints as morphological constraints. If a morphological constraint is not chosen properly, 
the model may converge to a physically implausible solution. 

On the other hand, we can sometimes make an assumption, not being certain it is 
correct. In case of degeneracy, a solution based on an assumption may be better than having 
no solution at all. One assumption might be the age of the coeval components. Assuming 
a particular type of evolutionary track, the luminosities from stellar evolution models may 
then be obtained (Pols et al. 1995). Another such assumption could be the distance to the 
binary, e.g. from astrometry. Yet another assumption may be that either or both components 
are main-sequence stars. Since a signihcant percentage of all stars are on the main-sequence, 
there is a fair chance that our assumption is correct. 

Applying main-sequence constraint to component (s) of the modelled binary means im¬ 
posing M-L-T-R relations for main-sequence stars (see e.g. Maikov (2003) for such relations 
specihc to EBs). Consequentially, given a single parameter (e.g. component’s effective tem¬ 
perature), all other parameters (its mass, luminosity and radius) are calculable. This in 
turn implies that, in case of circular and nearly-circular orbits, effective potential of the con¬ 
strained component is fully determined. Main-sequence constraint may be used for testing 
whether either or both stars may plausibly be main-sequence stars: depending on behavior 
of the value, such hypothesis may be accepted or rejected. 

Such additional constraints are not as straight-forward as was the case with morpholog¬ 
ical constraints. For example, by implying the condition: let the modeled binary be a main- 
sequence binary, we break the degeneracy by selecting the one solution that corresponds to 
that condition. This is why PHOEBE refers to these constraints as conditional constraints 
(CC). It is very important to emphasize that using conditional constraints improperly may 
lead to creating and propagating a circular argument: EBs provide absolute parameters for 
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stars, which can then be used to establish various calibrations. Conditional constraints on 
the other hand use calibrations to constrain derived parameters. Conditionally constrained 
solutions should thus never be used to establish calibrations of any kind. 

Recall from Fig. 6 that the solution from photometric and RV observations of our test 
binary indeed suffered from degeneracy in inclination and potentials. If we conditionally 
constrain both modelled components with the main sequence constraint, potentials and 
r22 are calculable and thus exactly known. The variation in their values is only a consequence 
of the variation in either of the main-sequence parameters (M, L, T or R) that accomodate 
for different orbital inclinations during the fit. Fig. 9 shows convergence tracers for a sim¬ 
ilar NMS heuristic scan as in Fig. 6, this time for f2i-r22 cross-section without (left) and 
with (right) the main-sequence constraint imposed on the model. Since the main-sequence 
constraint is very strong, there is no practical need for heuristic scan or parameter kick¬ 
ing (both hlj’s are calculable for the given inclination and convergence is thus assured from 
practically any point in the hyperspace); Fig. 9 (right) still depicts both heuristic scan and 
consecutive parameter kicks for comparison between convergence tracer shapes and slopes 
of unconstrained and main-sequence constrained model. It is evident that both solutions 
intersect, yielding the right solution. This is of course expected, since our test binary is in 
fact composed of two main-sequence components. 

One would hope that total eclipses reduce the degeneracy, but this also does not nec¬ 
essarily happen. If stars have comparable sizes (as is the case of our test binary), the 
duration of the eclipse totality is very short and limb darkening may obscure its flatness. 
Since geometrical parameters in case of total eclipses are better constrained (the correspond¬ 
ing hyperspace cross-sections feature very narrow valleys), parameter kicking may work to 
our disadvantage, knocking the solution far from the minimum by only a small parameter 
displacement. This issue will be addressed in detail in the follow-up paper. 

By using conditional constraining, we select a preferred subspace of model solutions. 
This is why extra care should be taken for the choice of adopted CCs. 


4.4. Interstellar and atmospheric extinction 

Although interstellar extinction has been discussed in many papers and quantitatively 
determined by dedicated missions (lUE, 2MASS, and others), the approach for EBs is often 
inadequate. Reddening is usually calculated by analytic approximation (e.g. Lang et al. 
1992) or extinction tables (e.g. Schlegel et al. 1998), using EB’s galactic coordinates and 
inferred distance to the binary. Such calculations are performed only for a single, effective 
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Fig. 9.— Convergence tracers for f2i-f22 cross-section without (left) and with (right) the 
main-sequence constraint imposed on the model. Similar to Fig. 6, panels (a) through (d) 
denote successive number of kicks (zero to three) and cross-hairs mark the correct solution. 
Comparing these results immediately shows that the intersection of the two attractors yields 
the correct solution. Note that there are two intersections because of the model symmetry 
to the labeling of the two components (primary and secondary roles of both stars may be 
interchanged). 



























wavelength of the given passband; the obtained value is then subtracted uniformly from all 
photometric observations in that passband. 

If the amount of reddening is not negligible and the binary components have signihcantly 
different surface temperatures, effective temperature of the binary is a function of phase 
because of eclipses. This is why in case of strongly reddened EBs it cannot be assumed that 
the correction due to reddening is well approximated by simply subtracting a constant. A 
dedicated study of this problem has been presented by Prsa & Zwitter (2004), here we only 
overview the conclusions. 

PHOEBE uses already described synthetic SED data for rigorous reddening corrections. 
It builds an intrinsic effective spectrum of the binary by Doppler-shifting and convolving 
the spectra of individual components as a function of phase. This intrinsic spectrum is 
then rigorously (wavelength-by-wavelength) reddened by the formula proposed by Cardelli 
et ah (1989) and convolved with the given passband transmission function. Photometric 
magnitude is then obtained by integrating the reddened spectrum over the given passband. 

There are two major implications of this improved scheme over the traditional constant 
subtraction, that are depicted in Fig. 10; 1) using effective passband wavelength to calculate 
the reddening constant introduces a signihcant systematic error in reproduced magnitudes. 
Effective wavelength of the passband is irrelevant for the reddening: it is the spectrum 
integral over the passband interval which must be the same in both approaches. 2) Even if 
the constant was calculated properly (by making sure that the integrals over the passband 
are the same), there would still be a measurable offset in both eclipses due to the change 
in effective temperature of the binary. This effect gains on significance as the temperature 
difference between both components grows and may reach ~ 0.2 mag or more in case of 
symbiotic binaries (Prsa & Zwitter 2004). 

Atmospheric extinction has a similar effect on photometric observations, since it also 
reddens the data as a function of wavelength. The study has shown that interstellar extinc¬ 
tion dominates the whole wavelength range in cases of strong reddening, while atmospheric 
extinction dominates the blue parts of the spectrum in cases of low-to-intermediate reddening 
(Prsa & Zwitter 2004). The main difference between interstellar and atmospheric extinction 
is that the latter is usually taken into account during initial reduction of the photometric 
data, so prior to running PHOEBE. 
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5. Conclusion 

This paper has presented current stage of an ongoing effort to unify proven ideas and new 
approaches of the EB field. PHOEBE is built on top of the WD code and decades of experience 
put into its development, trying to add its own pieces to the puzzle: new minimization 
algorithms, heuristic scans, parameter kicking and an extended set of physical constraints. 
It is continuously growing and maturing due to constructive feedback of many individuals. 
In future, PHOEBE aims to broaden its scope on all areas mentioned in this paper: numerical, 
scientific and technical. We conclude this paper by naming some of the goals PHOEBE has 
yet to achieve. 

1. Full-scale testing. PHOEBE’s scientific core is now ready for extensive testing on real 
data. Individual stars as well as large survey databases are ideal testing grounds to 
hunt down problems and improve those aspects that may now be lagging behind. 

2. Scripting. Although the graphical user interface (see Appendix A.2) is well suited for 
individual targets, its usability is very limited when the number of EBs is large. We 
must prepare for the upcoming missions like Gaia, since the shear number of observed 
EBs will be several orders of magnitude larger than the number of all already solved 
EBs of today or, for that matter, the number of astronomers in the world to solve 
them all, one-by-one. It is naive to believe that our procedures are already optimal 
and applicative to all sorts of EBs that are out there. 

3. New physics. Some of the ideas already mentioned in this paper are in their infancy. 
The SED must evolve into a consistent and reliable data source that enables us to 
not only confirm or reject the otherwise obtained solution, but to extract parameters 
from the spectra themselves. Once the SEDs are fully integrated in solution seeking, 
LD coefficients will have become obsolete, for intensities will then be computable from 
spectra and the LD effect will come out naturally. Individual components may be 
intrinsically variable and common types of variabilities may easily be recovered from 
the model (see Dallaporta et ah (2002) for an example of a (5-Sct companion). 

4. New numerical algorithms. By the ever-growing computer power, better and more 
powerful numerical algorithms are surfacing. Two very promising candidates are al¬ 
ready in testing: Adaptive simulated annealing (Ingber 1996) and Powell’s direction 
set method (Acton 1990). Both are based only on function evaluations, not numerical 
derivatives. 

5. New technical enhancements. With continuous help and support from users sharing 
their opinions and suggestions on PHOEBE discussion mailing lists, we are able to form a 
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wish-list and implement most needed features. Custom user-supplied passband trans¬ 
mission functions must be supported to enable the data obtained by any instrument 
and any hlter-set to be processed. 

Since PHOEBE is free (released under the GNU GPL - see Appendix), everyone with 
enthusiasm and interest may join in on this project! PHOEBE will keep improving. 


The authors would like to express their utmost gratitude to Robert E. Wilson, for 
spending numerous hours commenting and criticising the manuscript and for his continuous 
encouragement. We are also indebted to the referee of the paper, for making valuable 
suggestions that signihcantly improved the paper’s layout and clarity. Fruitful discussions 
with Dirk Terrell, Michael Bauer, Walter Van Hamme, Michael Sallman and Ulisse Munari 
throughout PHOEBE development are very much appreciated. Our thanks go to all PHOEBE 
users out there, supporting our work with constructive feedback. We would also like to thank 
Phillip J. Flower for his swift reply on updated coefficients given in Table 2. This work is 
supported by the Slovenian Ministry for High Education, Science and Technology. 


A. Appendix: technical information 

PHOEBE is released under the GNU General Public License (GPL) and is freely avail¬ 
able for download from http: //phoebe. f iz. uni-1 j . si. The package comes with thorough 
documentation: Reference manual. Tutorial and Application Programming Interface. For 
convenience, three discussion mailing lists are available to help users communicate and share 
opinions, ideas and enhancements to PHOEBE. The freedom of GPL enables anyone with 
interest to join in on future development. 


A.l. The back-end: PHOEBE scripter 

In its core, PHOEBE is a scripting language. This means that the user communicates with 
the program interactively by passing particular statements to perform particular actions. 

PHOEBE language is based on formal, context-free LALR(l) grammar. This means that 
strict and consistent grammar rules of scanning, parsing and evaluating user input are im¬ 
posed to achieve full support for arithmetics, nested loops, conditionals and function deh- 
nitions (see Aho et al. (1986) for specihc properties of LALR(l) grammar). It is written in 



ANSI C, which makes it portable to virtually any platform regardless of the operating system 
used. 

PHOEBE scripter consists of three layers. The lower-most layer is the WD code, the layer 
above is PHOEBE’s extension layer and the topmost layer is the interpreter. The underlying 
WD code is unchanged, which makes adaption to any future WD versions trivial. The 
extension layer contains scientihc add-ons that enhance basic WD applicativity. Finally, 
the interpreter’s purpose is to communicate with the user. Its plug-in awareness allows 
miscellaneous technical enhancements to be easily incorporated - the graphical user interface 
(GUI), main-sequence calculators etc. 


A.2. The front-end: PHOEBE graphical user interface 

All novel PHOEBE features discussed in the main part of the paper are implemented in 
PHOEBE’s back-end engine, the scripter. Although it is gratifying to achieve advancements 
in scientihc and numerical approaches, technical details that make scientist’s life easier are 
all-too-often overlooked. Based on current PHOEBE users’ feedback, the most prominent 
enhancement PHOEBE brings into the held is neither numerical nor scientihc, it is technical: 
a graphical user interface (GUI). No longer is it necessary to spend hours or even days 
learning the technicalities of a particular code; PHOEBE features a full-hedged, hexible and 
heavily structured GUI that brings the ease of clicking, observing and monitoring the process 
of solution seeking to the user. 

The GUI. Any implementation of a front-end is inherently system-dependent, and so is 
PHOEBE’s interface. The GUI is designed to run under any Linux (or other Linux-compatible) 
operating system. It should be noted nevertheless that the GUI is merely a plug-in to PHOEBE 
scripter, so when a need for a diherent GUI on a diherent operating system arises, it is only 
a matter of building a front-end - the back-end will remain the same, portable to all ANSI-G 
compliant architectures. 

PHOEBE’s GUI is written with GTK+ graphical library, a free standard component of 
virtually any Linux of today. It consists of the main screen, the snapshot of which is depicted 
in Fig. 11. The main window is used for basic user interaction - changing parameter values, 
obtaining statistics on observations, plotting star hgures etc. From the button menu on 
the bottom of the main window users may open auxiliary windows. They are used to plot 
photometric light curves, RV curves, to initiate the fit or to write PHOEBE scripts. The 
interface is consistent with the rest of the operating system, so users with elementary Linux 



experience should have PHOEBE up and running in no time. 


Plotting observations and solutions, star figures. Plotting light curves and RV curves 
throughout the fitting process proves to be extremely useful for consistency checking. PHOEBE 
supports data binning, overplotting, phase aliasing and plotting 0 — C residuals. In addition, 
PHOEBE plots star figures as they appear on the sky (see the lower right panel of Fig. 1 for 
an example) in any given phase. It is also capable of enhancing the location of spots, making 
animated cartoons or e.g. depicting apsidal motion of eccentric binaries. 


LD interpolation functions. Another very important aspect of modeling EBs that was 
the sole responsibility of WD users are the limb darkening (LD) coefficients. Rather than 
trusting WD to retrieve the values of the coefficients by numerical fits, it is much safer to 
use precomputed LD tables, e.g. Van Hamme (1993). PHOEBE can retrieve correct values of 
these coefficients from external LD tables at each iteration step. This speeds up convergence, 
improves consistency of the derived solution and avoids the need for manual adjustment 
during the iteration process. 


Measuring parameter correlation, handling degeneracies. Throughout the paper 
we have discussed the dark world of parameter correlations and degeneracies. It is crucial 
for any user to be aware of these problems. An excellent review on what one may expect from 
light curve modeling is given by Wilson (1994), which is both informative and entertaining, 
arming the user against the caveats that await him. This is why in addition to WD’s 
correlation matrix (Wilson & Van Hamme 2003) PHOEBE plots parameter histograms and 
convergence tracers during the fit. If we are able to see and inspect the solution in each 
iteration step and, perhaps more importantly, if we are able to play with the solution as we 
see fit, we will have more fun with it and more consistent physics is bound to emerge. 
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Fig. 10.— Johnson B magnitnde difference between reddened and nnreddened observations. 
£^(B — V) = 1 and i = 90° is assnmed with other parameters listed in Table 1. Left: The 
discrepancy between the rigoronsly applied reddening (points) and the constant snbtraction 
approach (solid line). The difference dne to erroneons approach is ~ 0.06 mag. Snbtracted 
constant was obtained from the effective wavelength (Agg = 4410.8 A) of the Johnson B 
passband transmission curve. Right: Overplotted light curves with the subtraction constant 
calculated properly in out-of-eclipse regions. There is still a measurable difference of ~ 
0.01 mag in eclipse depth of both light curves. Adapted from Prsa & Zwitter (2004). 











33 



Fig. 11.— A snapshot of PHOEBE graphical user interface in action. 































